Variable Koeffizienten¶
Das folgende Beispiel folgt dem analogen aus dem interaktiven Kurs von Joachim Schöberl [SHL21].\(\DeclareMathOperator{\opdiv}{div}\) \(\DeclareMathOperator{\setR}{R}\)
from ngsolve import *
from ngsolve.webgui import Draw
Problem¶
Ein Setup mit unterschiedlicher Wärmeleitungskoeffizienten wir mit Hilfe der Gleichung
modelliert, wobei \(\lambda=\lambda(x)\) der ortsabhängige Wärmeleitungskoeffizient ist. Der zugehörige Wärmefluss
ist gegeben durch den Gradient der Temperatur \(\nabla u\). Im Fall, dass \(\lambda\) diskontinuierlich ist, wird die Gleichung im Sinne von Distributionen verstanden. Dies beinhaltet Interface Bedingungen: Die Temperatur auf der linken und rechten Seite sind gleich und der Wärmefluss der linken in die rechte Seite müssen gleich gross sein:
Die variationelle Form des Problems ist gegeben durch: finde \(u \in H^1(\Omega)\) mit
Diskontinuierliche Koeffizienten bilden kein Problem. Beide Interface Bedingungen werden erfüllt:
Stetigkeit der Temperatur \(u\) durch die Stetigkeit der Ansatzfunktionenraumes
Stetigkeit des Wärmeflusses in einem schwachen Sinne, durch Neumann Randbedingungen.
Geometrie und Mesh¶
from netgen.geom2d import *
geo = SplineGeometry()
geo.AddRectangle( (0,0), (1,1), leftdomain=1, rightdomain=0,
bcs=['b','r','t','l'])
geo.AddCircle( (0.3,0.7), 0.1, leftdomain=2, rightdomain=1)
geo.AddRectangle ( (0.2,0.2), (0.9,0.3), leftdomain=3, rightdomain=1)
geo.SetMaterial(1, "air")
geo.SetMaterial(2, "source")
geo.SetMaterial(3, "bar")
mesh = Mesh(geo.GenerateMesh(maxh=0.03))
mesh.Curve(3)
Draw (mesh,)
BaseWebGuiScene
Die Subdomains (Materials) sind gegeben durch:
print (mesh.GetMaterials())
('air', 'source', 'bar')
Die Geometrie hat folgende Ränder (innere und äussere):
print (mesh.GetBoundaries())
('b', 'r', 't', 'l', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default')
V = H1(mesh, order=3, dirichlet="b|r")
u = V.TrialFunction()
v = V.TestFunction()
Wir definieren nun das Material in dem wir den Wärmekoeffizient stückweise konstant ansetzen:
lamvalues = { "air" : 1, "bar" : 1e2, "source" : 2 }
lam = CoefficientFunction(
[lamvalues[mat] for mat in mesh.GetMaterials()])
Draw (log(lam), mesh, "log lambda")
BaseWebGuiScene
Wir nützen die Wärmeleitfähigkeit für die Definition der Bilinearform:
a = BilinearForm(V)
a += lam*grad(u)*grad(v)*dx
f = LinearForm(V)
f += 1*v*dx("source")
a.Assemble()
f.Assemble()
<ngsolve.comp.LinearForm at 0x7fb5d7e2d430>
und lösen das System:
gfu = GridFunction(V)
gfu.vec.data = a.mat.Inverse(V.FreeDofs()) * f.vec
Für die Wärmeverteilung erhalten wir:
Draw (gfu, mesh, "temperature")
BaseWebGuiScene
Für den Gradient der Wärmeverteilung \(\nabla u\)
Draw (grad(gfu), mesh, "gradient",vectors=grad(gfu))
BaseWebGuiScene
und für den Wärmefluss \(-\lambda \nabla u\):
Draw (-lam*grad(gfu), mesh, "heatflux")
BaseWebGuiScene